﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.FFTPack
{
    /// <summary>
    /// Sine FFT transform of a real odd sequence. 
    /// </summary>
    [Serializable]
    public class RealDoubleFFTOdd : RealDoubleFFTMixed
    {
        /// <summary>
        /// This field holds the norm factor.
        /// </summary>
        private double normFactor;

        /// <summary>
        /// This field holds the wave table.
        /// </summary>
        private readonly double[] wavetable;

        /// <summary>
        /// This field holds the dimension.
        /// </summary>
        private readonly int ndim;

        /// <summary>
        /// Initializes a new instance of the <see cref="RealDoubleFFTOdd"/> class.
        /// </summary>
        /// <param name="n">The size of a real data sequence. When (n+1) is a multiplication of small numbers 
        /// (4, 2, 3, 5), this FFT transform is very efficient.</param>
        public RealDoubleFFTOdd(int n)
        {
            ndim = n;
            normFactor = 2 * (n + 1);
            int wtable_length = 2 * ndim + ndim / 2 + 3 + 15;
            if (wavetable == null || wavetable.Length != wtable_length)
            {
                wavetable = new double[wtable_length];
            }
            Sinti(ndim, wavetable);
        }

        /// <summary>
        /// This field can be used to normalize this FFT transform. This is because a call of forward 
        /// transform followed by a call of backward transform will multiply the input sequence by the 
        /// norm factor.
        /// </summary>
        /// <value>The norm factor to use for the transform.</value>
        public double NormFactor
        {
            get { return normFactor; }
            set { normFactor = value; }
        }

        /// <summary>
        /// Forward sine FFT transform. It computes the discrete sine transform of an odd sequence.
        /// </summary>
        /// <param name="x">An array which contains the sequence to be transformed. After FFT, 
        /// x contains the transform coeffients.</param>
        public void ForwardTransform(double[] x)
        {
            Sint(ndim, x, wavetable);
        }

        /// <summary>
        /// Backward sine FFT transform. It is the unnormalized inverse transform of ForwardTransform();
        /// </summary>
        /// <param name="x">An array which contains the sequence to be transformed. After FFT, x contains 
        /// the transform coeffients.</param>
        public void BackTransform(double[] x)
        {
            Sint(ndim, x, wavetable);
        }

        private void Sint1(int n, double[] war, double[] wtable)
        {
            double sqrt3 = 1.73205080756888;
            int modn, i, k;
            double xhold, t1, t2;
            int kc, np1, ns2;
            int iw1, iw2, iw3;
            double[] wtable_p1 = new double[2 * (n + 1) + 15];
            iw1 = n / 2;
            iw2 = iw1 + n + 1;
            iw3 = iw2 + n + 1;

            double[] x = new double[n + 1];

            for (i = 0; i < n; i++)
            {
                wtable[i + iw1] = war[i];
                war[i] = wtable[i + iw2];
            }
            if (n < 2)
            {
                wtable[0 + iw1] += wtable[0 + iw1];
            }
            else if (n == 2)
            {
                xhold = sqrt3 * (wtable[0 + iw1] + wtable[1 + iw1]);
                wtable[1 + iw1] = sqrt3 * (wtable[0 + iw1] - wtable[1 + iw1]);
                wtable[0 + iw1] = xhold;
            }
            else
            {
                np1 = n + 1;
                ns2 = n / 2;
                wtable[0 + iw2] = 0;
                for (k = 0; k < ns2; k++)
                {
                    kc = n - k - 1;
                    t1 = wtable[k + iw1] - wtable[kc + iw1];
                    t2 = wtable[k] * (wtable[k + iw1] + wtable[kc + iw1]);
                    wtable[k + 1 + iw2] = t1 + t2;
                    wtable[kc + 1 + iw2] = t2 - t1;
                }
                modn = n % 2;
                if (modn != 0)
                {
                    wtable[ns2 + 1 + iw2] = 4 * wtable[ns2 + iw1];
                }
                Array.Copy(wtable, iw1, wtable_p1, 0, n + 1);
                Array.Copy(war, 0, wtable_p1, n + 1, n);
                Array.Copy(wtable, iw3, wtable_p1, 2 * (n + 1), 15);
                Array.Copy(wtable, iw2, x, 0, n + 1);
                Rfftf1(np1, x, wtable_p1, 0);
                Array.Copy(x, 0, wtable, iw2, n + 1);
                wtable[0 + iw1] = 0.5 * wtable[0 + iw2];
                for (i = 2; i < n; i += 2)
                {
                    wtable[i - 1 + iw1] = -wtable[i + iw2];
                    wtable[i + iw1] = wtable[i - 2 + iw1] + wtable[i - 1 + iw2];
                }
                if (modn == 0)
                {
                    wtable[n - 1 + iw1] = -wtable[n + iw2];
                }
            }
            for (i = 0; i < n; i++)
            {
                wtable[i + iw2] = war[i];
                war[i] = wtable[i + iw1];
            }
        }

        private void Sint(int n, double[] x, double[] wtable)
        {
            Sint1(n, x, wtable);
        }

        private void Sinti(int n, double[] wtable)
        {
            double pi = Math.PI;
            int k, ns2;
            double dt;

            if (n <= 1)
            {
                return;
            }
            ns2 = n / 2;
            dt = pi / (double) (n + 1);
            for (k = 0; k < ns2; k++)
            {
                wtable[k] = 2 * Math.Sin((k + 1) * dt);
            }
            Rffti1(n + 1, wtable, ns2);
        }
    }
}